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Abstract 

We present a new empirical pseudopotential (EPM) calculation approach to simulate the million 
atom nanostructured semiconductor devices under potential bias using the periodic boundary con- 
ditions. To treat the non-equilibrium condition, instead of directly calculating the scattering states 
from the source and drain, we calculate the stationary states by the linear combination of bulk 
band method and then decompose the stationary wave function into source and drain injecting 
scattering states according to an approximated top of the barrier splitting (TBS) scheme based on 
physical insight of ballistic and tunneling transport. The decomposed electronic scattering states 
are then occupied according to the source/drain Fermi-Levels to yield the occupied electron den- 
sity which is then used to solve the potential, forming a self-consistent loop. The TBS is tested 
in an one-dimensional effective mass model by comparing with the direct scattering state calcula- 
tion results. It is also tested in a three-dimensional 22 nm double gate ultra-thin-body field-effect 
transistor study, by comparing the TBS-EPM result with the non-equilibrium Green's function 
tight-binding result. We expected the TBS scheme will work whnever the potential in the barrier 
region is smoother than the wave function oscillations and if it does not have local minimum, thus 
there is no multiple scattering as in a resonant tunneling diode, and when a three-dimensional 
problem can be represented as a quasi-one-dimensional problem, e.g., in a variable separation ap- 
proximation. Using our approach, a million atom non-equilibrium nanostructure device can be 
simulated with EPM on a single processor computer. 

PACS numbers: 
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I. INTRODUCTION 



According to the roadmap of the Semiconductor Industry Association 1 , MOSFET (metal 
oxide semiconductor field-effect transistor) channel length will scale down to 22 nm by 
2012. In such nanosized devices, quantum mechanical effects play a big role in determining 
the properties of the system. New quantum mechanical features, like the fact that the 
electron mean free path is larger than the device dimensions and the single quantum state 
levels, can be used to enhance device performance and form new functionalities^ - -. On 
the other hand, as the size reduces, new obstacles emerge^, such as the short channel 
effects, source/drain off-state quantum tunnelling current, barrier current leakage and single 
dopant random fluctuation^. Over the past 20 years, many methods have been developed to 
incorporate the quantum mechanical effects into the device simulation^ - -. The first of that 
is the inclusion of some quantum mechanical effective potentials in the drift equation based 
on the gradients of the charge density^* 1 ^. Such gradient terms make the charge density 
smooth near the Si/Si02 interface, hence incorporating some of the quantum mechanical 
effects. The so called quantum Poisson drift equation, or quantum hydrodynamic model have 
been extensively used for device simulations. However, as the device size shrunk further, it 
was realized that the quantum mechanical wave functions need to be calculated explicitly. 

There are many ways to do the Poisson-Schrodinger's equation depending on the problem 
to be studied and the computational costs. One way is to calculate the quantum mechani- 
cal local density of states, then apply Boltzmann transport equation based on such density 
of states 1 ^ 1 ^. However, in the ballistic size region, the use of Boltzmann equation itself is 
questionable. In such cases, the direct solution of the open boundary condition scattering 
states based on the Schrodinger's equation is necessary. For example, this has been done 
for the ID cases using the nonequilibrium Green's function approach base on tight-binding 
model 14 . The use of Green's function also provides a way to incorporate the inelastic scat- 
tering processes in the formalism. Recently, three dimensional devices models with hundreds 
of thousands of atoms have been simulated using the tight-binding model based on the cal- 
culation of scattering states 1 ^^. But thousands of computer processors are needed for such 
direct 3D simulations. There are also effective mass calculation for 2D systems using the 
scattering state approach™^. It involves the solutions of linear equations in the dimension of 
the number of real space grid points. Overall, the direct simulation for the 3D device model 
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based on quantum mechanical transport equation remains to be nontrivial. Thus, it will be 
very useful if there is a faster way to do the simulations. One possible approach is to use 
the stationary eigen states of a closed system (e.g., periodic system) Schrodinger's equation 
to represent the quantum mechanical effects of the scattering states^ - — . The calculation of 
the eigen states for a closed boundary condition (e.g., periodic boundary condition) prob- 
lem is much faster than the calculation of an open boundary condition scattering states. 
This approach is plausible in the sense, at the zero bias potential between the source and 
drain, the open boundary condition problem is the same as the closed boundary condition 
problem. Thus, the closed boundary condition solution is a good starting point for the open 
boundary condition problem. Many of the quantum mechanical effects have already been 
represented at this close boundary condition level. The challenge however is to find a good 
approximation to get the charge density from the stationary wave functions when there is 
a large bias potential, and to find the corresponding electron current. The possibility of 
such an approximation relies on the fact that the potential profiles for many semiconductor 
electronic devices are often relatively simple and smooth, when there is no local minimum 
in the potential, thus no pseudo localized states, the perspective of the coherent multiple 
scattering is small. Thus, we will exclude ourselves from cases like the resonant tunneling 
diode (RTD), or complicated molecular electronics. The approximation might be feasible 
especially in the regime of ballistic transport, here we define it not only as elastic transport 
(ignoring electron-phonon scattering), but also as a current overcomes a barrier, then flushes 
through down hill without multiple scattering^ 1 ^. In the down hill flushing regime, sim- 
ple approximations like the WGK approximation can be applied. Note that, the scattering 
states satisfy the same Schrodinger's equation as the stationary states, albeit their boundary 
conditions are different. Thus the eigen states might contain the information as needed for 
a transport problem (e.g., the local density of states). In this work, we will present a new 
way to obtain the scattering states from the stationary eigen states. As can be seen in the 
following sections, our approach is physically intuitive, easy to implement, and tested to be 
accurate for the smooth barrier potential cases as described above. 

Besides the open boundary condition and close boundary condition problem, there is an 
issue of what Hamiltonian to be used to describe the electronic property of the system. We 
will use empirical pseudopotential method (EPM) as our Hamiltonian to study the problem. 
We have mentioned the tight-binding modeU^&2£>2& and effective mass modet£&^ above. 
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Another often used model is the the k ■ p model which can be used to describe the multiple 
valence band state s 18126 ' 27 '. However, as shown by Esseni and Palestri et al.— , the k ■ p can 
significantly misrepresent the electron density of states (DOS) for a Si inversion layer, and 
the indirect band gap nature of bulk Si presents real challenges for the k ■ p/effective-mass- 
like models 2 ^. It has been shown that in some cases, the use of the full band structures is 
important in simulating Si nanodevices 2 ^ 3 ^. The EPM is an accurate method to describe 
the full semiconductor band structures and electron wave functions. Within EPM, the total 
electron potential V(r) is described as a superposition of the spherical atomic EPM potentials 
v a (|r|) as V(r) = ^2 R v a (\r — R\), while f a (|r|) for atom type a is fitted to experimental bulk 
band structures, and R is the atomic positions. The wave functions in the EPM approach 
are expanded using plane wave basis sets. The atomic feature of EPM is important for 
simulating small nanodevices where the single atomic characters become significant. EPM 
can also describe other effects like strain, heterostructure, and semiconductor alloying and 
components, which are the current research topics for Ge/Si and InAs devices. Another 
reason to use EPM is one particularly fast algorithm to calculate the electron eigen states 
in a periodic boundary condition situation. This is the Linear combination of bulk band 
(LCBB) method 31 . The LCBB uses the bulk band states instead of the original plane 
waves as the basis set to expand the electron wave functions. As a result, the number of 
the basis function can be truncated to be less than 10,000 by selecting a finite number 
of k-points and band index. The resulting Hamitlonian matrix can be diagonalized by a 
single processor computer within an hour. The LCBB eigen energy is within 10 meV of 
the directly calculated results using the plane wave basis set. The LCBB is an atomistic 
method since each individual atom can be replaced by another atom, and it can describe 
the strain effects 3 ^. The computational time of LCBB method is roughly independent of 
the system size (since it depends mostly on the number of Bloch basis functions which in 
many case is independent of the system size), thus can be used to calculate million atom 
nanostructures^ 2 -. 

II. SIMULATION APPROACH 

The simulation approach follows the typical iterative scheme which solves the Schrodinger 
equation (1) and Poisson equation (2) self-consistently. 
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(1) 



V[e(r)V<£(r)] 



47i[p(r) -n(r) + Nj 



(r)-N-(r)], 



(2) 



here V(r) = ^2jtV a {r — R) is the total empirical pseudo-potential of silicon crystal 3 ^. 
V str (r) is the quantum confinement potential representing the Si02 barrier layer and the 
potential well between the source and drain for the artificial periodic boundary condition. 
<p(r) is the electrostatic potential which is solved by the Poisson equation (2). e(r) is the 
position dependent dielectric constant; N£(r) and N~(r) are donor and acceptor nuclei 
charges, to be treated as continuous charge densities in the current calculation; p(r) is 
the hole charge density while n(r) is the occupied electron charge density. A more detail 
description of the device setup and the way to solve the Poisson equation have been described 



As mentioned before, the main issue to be addressed in the current paper is to calculate 
occupied electron density n(r) and the total current from eigen state pairs {Ei,ipi(r)} of 
Eq.(l). Previously 33 ' 34 , we have used the WKB approximation as a weight function as well 
as the local Quasi-Fermi Potential to solve this problem. In the following subsections we 
will present a new approach based on physical intuition, and to be tested numerically. It is 
also stable in both ballistic and tunneling cases. Compared to our previous approaches 3 ^ 4 -, 
the present approach gives a smoother charge density and current^ 3 -, and is conceptional 
in the ballistic regime instead of the thermal scattering regime^ 4 -. More importantly, our 
numerical tests show that the results are surprisingly accurate compared to open boundary 
condition scattering state calculations. In this paper, our formalism will be presented in 
a heuristic style, instead of rigorous derivations from the original scattering state problem. 
They are derived based on a few simple principles and assumptions. They are tested for 
typical systems representing the problems we intend to solve. As will be discussed below, our 
top of the barrier splitting (TBS) algorithm will be based on an one-dimensional effective 
mass algorithm, while our original eigen-state wave function could be calculated by EMP 
(e.g., in subsection B). Thus it represents a hybrid approach. The effective mass TBS will 
not devalue the final result to the effective mass level since the features of the atomistic EPM 
calculation, e.g., the atomistic wave function, the non-parabolicity of the kinetic energy, the 



in Ref 34. 
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multiple valley, will be retained in the overall procedure. 

Our problem can be summarized as how to use {Ei, ipi{r)} to occupy the system with a 
left and a right Fermi energies Ep and Ep to get the occupied charge density n(r) and how 
to estimate the current /. Formally, this problem can be solved by occupying the scattering 
states ipg(r, E),ips(r, E) by the left and right Fermi energies respectively, 

n(r) = j (|V£(r, E)ff(E - E L F ) + \^(r, E)\*f(E - ^))(^)~ 1 dB (3) 

I = J (T s L (E)f(E - E L F ) - T s R (E)f(E - Ep))(^^)~ 1 dE (4) 

where f(x) is the Fermi-Dirac distribution function for a given temperature T, and T^(E), 
T^(E) are the left and right scattering states transmission coefficients. The scattering states 
are wave functions satisfying the same Schrodinger's equation as in Eq.(l), but with an open 
boundary condition 40 . So, our question is whether we can use {Ei, ipi{r)} to mimic the effects 
of the scattering states ipg(r, E) and i^§(r, E), or at least their energy integrated properties 
n(r) and /. 

Subsection A will describe the one-dimensional splitting algorithm while subsection B 
will extend it to three-dimensional case which will then be incorporated with the LCBB 
calculation for our device simulation. 



A. One-Dimensional Model 



Suppose that the electron is running along the x direction, here we will base our formalism 
on an ID effective mass Hamiltonian: 

{~^ + V(x))Ux) = m{x) (5) 

We will assume V(x) is a smooth function as shown in Fig.l, and it has a barrier at the 
center just like the situation in a transistor, more specifically, the variation is slower than 
the wave function oscillations, and there is no local potential minimum in the barrier region. 
Now, for each eigen wave function ipi{ x ) under a closed boundary condition (in our case, a 
periodic boundary condition), we will like to break it into the left injecting ij)f(x) and the 
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right injecting ip^(r) scattering states. We first require the scattering states satisfying the 
charge conservation rule: 



\^{x)\ 2 = \tf{x)\ 2 +\^{x)\ 2 



(6) 



The above requirement Eq.(6) is needed so that at equilibrium (Ep = Ep), occupying 
\^i(x)\ 2 and \ipf-(x)\ 2 is the same as occupying \il)i(x)\ 2 as in the original closed system 
problem. The occupied electron density as well as the current density will be obtained from 
equation (3) and (4). 

Now for a ID potential V(x) shown in Fig.l, we will have a unique maximum barrier 
height V m at x m (since there is no potential minimum), and left source potential Vl and 
right drain potential Vr as shown in Fig.l. In the following we will distinguish three cases: 
ballistic case with Ei > V m , tunneling case with Vl < Ei < V m , and stationary case of 
V R <Ei< V L . 

For the ballistic case, we require that i^f{x) and ip^(x) are the same at x m , i.e., 



This requirement is necessary, as we will see below, coupled with the current model, 
this will guarantee the current of the left and right scattering states equal. Such equality 
is needed so at equilibrium when both \ipf(x)\ 2 and {^^{x)^ are occupied, there is no net 
current. 



\^i{x)\ 2 , the current becomes ballistic (flushing down hill) for x > x m . Now we will use 
a ballistic approximation, where the current equals Vi(x)\ipf(x)\ 2 . Since the current must 
be a constant, independent of x, it thus must be equal to Vi{x m )\^f (x m )\ 2 . Similarly, for 
the right injecting wave \ip^(x)\ 2 , the current becomes Vi(x)\ip^(x)\ 2 = Vi(x m )\ipf(x m )\ 2 for 
x < x m . This leads to: 



\tf(x m )\ 2 = \lfM\ 2 = \\Mx m )\ 2 



(7) 



We define a transport velocity as Vi(x) = 




1 ; Wl . Now for the left injecting wave 

m x 





(8) 



and 
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mx)\ 2 = 



\A{x)\ 2 - V ^\U^m)\ 2 (X>x m ) 

(9) 



Note, the first equation in Eq. (8) and second equation in Eq. (9) comes from the current 
conservation law and the ballistic expression for the current, while the second euqation in 
Eq.(8) and first equation in Eq.(9) come from Eq.(6). Now, the currents for the left and 
right scattering states equal to: = Jf 1 = \vi(x m )\ipi(x m )\ 2 . Note that, the first equation 
in Eq.(8) and second equation in Eq.(9) can also be derived from WKB approximation. 

For the tunneling case, the eigen energy E { cross the barrier potential V(x) at x L and 
xr, xl < x m < x R as shown in Fig.l. 

Now for the region xl ^ x ^ x R , we consider the barrier interval (x,x' m ) (same for 
(x' m ,x)), here x' m is the minimum position for |^>j(a;)| 2 within the retion [xl,xr] (note that 
x m and x' m might not be the same, and x' m depends on the state index i). We can assume 
that the ratio of decay of \ip^(x)\ 2 from x to x' m is the same as the ratio of decay of \ip^(x) | 2 
from x' m to x. We also assume that the amplitude splitting equation of Eq.(7) holds at x' m . 
We thus have: 



m*)\ 2 \mx'j\ 



kvji 2 ~ mxw 



(10) 



i#(^ji 2 = i^vjr = i|^(x'ji 2 (ii) 

From the above two equations (10), (11), and Eq.(6), we can solve the i/jf' R (x) within 
the region as, 

\ti' R (x)\ 2 = lmx)\ 2 ± vm*)\ A -Mx'j\ A ) (12) 

where "+" is for the left injecting wave at x < x' m and the right injecting wave at x > x' m , 
while "-" is for the left injecting wave at x > x' m and the right injecting wave at x < x' m . 

Outside of the barrier (after the left and right injecting wave function tunnel out of the 
barrier), it is assumed that the current becomes ballistic. In that case, the wave function 
amplitude equals J/vi(x), here J is the tunneling current. All these lead to: 
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\^{ X )\ 2 = JK/Viix) 

|^(x)| 2 = |^(x)| 2 -|^(x)| 2 



(x < X L ) 
(x < X L ) 



(13) 



and 



\ri{x)\* = Jt/v t {x) 
|^(x)| 2 = |^(x)| 2 -|^(x)| 2 



(x > X R ) 
(x > X R ) 



(14) 



Note, the tunneling currents Jf, for the right and left scattering states should 
be the same (e.g., when they are both occupied, there will be no net current). We 
will use Jf- = J[ = ^\t/jf t (xL)\\ipi(xR)\^/2m*(V(x m ) — Ej). Another possible choice is: 
J? — Jt — \rnaxi\ijjf- (x L )\ 2 , \^(x R )\ 2 )^/2m* x (V(x m ) - E^). Although, these two choices 
look a bit arbitrary, in reality, we found no practical difference in our test between these two 
choices, we will thus use the former one. 

There seems a singularity problem at the classical turning points xl and xr in the first 
equation of Eq.(13) and Eq.(14). However, the charge density measure of this singular- 



ity is zero. More explicitly, Q_ t \^(x)\ 2 dx = J?yJ ^g^y f^_ e ^=dx = J t R ^ 



remove this singularity. Nevertheless, a small kink can be observed in Fig. 3. 

For the third case Vr < Ei < Vl, it is stationary, thus i>f{x) = ipi(x), ipi(x) = 0, and 



Now the occupied electron density as well as the total current can be evaluated as 



Once we have the electron eigen states {E^ip^x)}, the occupied electron density as well 
as the current / can be evaluated from the above model equations. We will call this the top 
of barrier splitting (TBS) model. 




J* = Jt = 0. 



n(x) = Ei \mx)?m - E L F ) + iVf (s)| 2 /(i* - E§) 
Jtot = Ei JifiEi - E L F ) - J?f(E t - E*) 



(15) 
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To test the validity of the above ID splitting model, we have calculated the scattering state 
exactly using the transfer matrix (TM) method 42 . For ID and effective mass Hamiltonian, 
this is easy to do. There is no numerical stability problem caused by the multiple band 
situation. We have chosen a test potential V(x) with a form as: 

V(x) = V L ( X < X\) 

< V(x) = -V g sin{^^) + (Vh - V L )£% + V L {x x <x< x 2 ) 

V{x) = V R {x> x 2 ) 

V 

(16) 

here X\ = 20 nm, x 2 = 45 nm, Vl — 0, Vr — —0.5 eV, and effective mass m* = 0.19 as 
for silicon. Note the shape of the potential can be modified by changing V g and Vl — Vr. 
We have tested other potential shapes (but with no potential minimum), and find similar 
results as shown below in terms of the accuracy of the TBS results. 

In Fig. 2 (a), the ballistic wave function splitting results are presented for two periodic 
boundary condition (by connecting x = nm potential to x = 65 nm potential) eigen states 
with eigen energies E 1 = 0.618 eV and E 2 = 0.620 eV. With V g = -0.8 eV, the V rn is at 
0.56 eV. Thus, both of these two states belong to the ballistic situation. Their decomposed 
\ip^(x)\ 2 are shown in Fig. 2(a) (the \il)f-(x)\ 2 looks similar from the opposite side). One 
can see that the constructed individual \ip^(x)\ 2 using Eq(6) can be negative. This might 
sounds alarming. However, the summation of these two scattering states (with very close 
eigen energies) |T/>f(o;)| 2 + [^(x)! 2 are all positive as shown in Fig.2(b). Furthermore, the 
summed result resembles closely to the TM calculated scattering state wave function at 
E = (Ei + E 2 )/2. Thus, as claimed previously, it is the energy integrated properties which 
resemble the directly calculated results, not the individual TBS scattering states. Note 
that, for ballistic case wave functions ^(x)! 2 with Ei > V m , the eigen states always come 
in pairs with very close eigen energies (rather like the sin and cos case). There is however 
and issue for low temperature occupation. What happens if | 2 is occupied while |V^| 2 
is not. Theoretically, this can be overcome by increasing the size of the source and drain 
regions, as a result, the difference between E\ and E 2 will diminish. In practice, at room 
temperature, we didn't find this is a problem. Tests for other ballistic case eigen states show 
similar behavior as in Fig. 2. 
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In Fig. 3, the decomposed (x)| 2 and \il)f(x)\ 2 scattering states for a tunneling case 
\ipi(x) \ 2 with E = 0.1814 eV are shown. As one can see, the decomposed state resemble 
closely the directly TM calculated scattering states. At x L for \i>f{x)\ 2 and x R for \^[(x)\ 2 , 
there is a kink to join the wave functions from Eq(12) and Eq(13) as discussed above. But 
in practice, since it happens in such a small amplitude, and the measure of this singularity is 
zero, it rarely matters. Other tunneling states show similar behavior. Note that for tunneling 
case, there is no pairing for the periodic boundary condition eigen states {Ei,ipi(x)}. 

In Fig. 4, the occupied charge densities calculated from Eq(15) are shown using Ep = 0.6 
eV, Ep = 0.1 eV, and V g = —0.8V at room temperature. We see that the TBS gives almost 
an identical occupied charge density as the TM directly calculated scattering state results. 
Note, here we are not doing any selfconsistent calculation yet, but such charge density is the 
first step towards the selfconsistent calculation (as will be done in the 3D calculation later 
in the paper). The calculated current for this case is 4.44 x 10~ 5 a.u. for the TBS model, 
while it is 4.58 x 10~ 5 a.u. for the TM directly calculated scattering state model, differs by 
only 3%. 

Now, we change V g and calculate the current I from left source to the right drain. This 
will give us the I-V curve. The results are shown in Fig.5. We see that the TBS model 
and the direct scattering state calculation yield almost the same results over a wide range 
of current amplitudes. This shows the accuracy of our TBS model for device simulations at 
least for one dimensional system. Note that, using Fermi-Diract distribution f(x) in Eq.(15), 
the small current in the region of V g < —0.9 eV comes mostly from the over-the-hill-top 
ballistic current caused by the ~ exp(-x) Boltzmann distribution. This is because as the 
barrier height increase, the tunneling current (below the hill top) decreases faster than the 
Boltzmann distribution exp(-x) for the over-the-hill-top states {ipi(x),Ei} to be occupied. 
To test the approximation on the pure tunneling current, we have also used an artificial 
occupation function f t (x) to replace the Fermi-Diract distribution in Eq(15). f t (x) = 1 for 
x < 0, ft(x) = cos 2 (xix/&) for < x < 4, and f t (x) = for x > 4. Thus, for large barrier 
height, there will be no over-the-hill-top ballistic current, instead all the current will come 
from the tunneling. From Fig.5, we see that this pure tunneling current also agree well 
between the TBS result and the direct scattering state calculation results. 
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B. Three-Dimensional Model 



In above, we have obtained an excellent ID algorithm. To extend that to three- 
dimensional case, we will first still base on an effective mass Hamiltonian, and employ 
the variable separation approximation between the x direction and y, z directions. Such 
separation is a good approximation in cases where there is a fast variation of the potential 
V(x,y,z) in the y, z directions, but a relatively slow variation in the x direction. This is 
applicable to many device systems with a narrow current channel. Such approximations 
have often been used to solve the 3D Hamiltonian eigen state problem. Here, we will use 
it to construct our splitting algorithm. Under the variable separation approximation, the 
three-dimensional wave function can be written as ip(r) = £(x)9(x,y,z) and 6{r) satisfies 
/ \9(r)\ 2 dydz = 1. We first require 

(-2^1? - + V '<*' »■ * Z > = * Z > (17 > 

This is a 2D eigen value problem with x as a parameter. Now, if we assume 9(x,y,z) 
varies slowly with x, thus its first and second order derivatives in x direction can be ignored, 
then satisfying the 3D Schrodinger's equation: 



Id 2 Id 2 Id 2 

( 77^7 7^7 7^7 + V(x, y, z))((x)9(x, y, z) = e{(x)8(x, y, z) (18) 

v 2m* dx 2 2m* y dy 2 2m* z dz 2 V ' y ' " SK J V ' y ' ; ; v >y> ; \ j 

is equivalent to satisfying the ID equation: 

{ -^ x & + E{xmx) = eax) (19) 
Now if by using other eigen state solvers (e.g., LCBB) we have obtained a three- 
dimensional eigen state pair {e, tp(r)} (we have droped the band index i for simplificty) , 
we can then write ((x) as: 

C 2 (x) = J mx,y,z)\ 2 dydz (20) 

and 9(r) as 9 2 (r) = ifj 2 (r)/( 2 (x). Once the one-dimensional effective potential E(x) is 
known, we can use the ID algorithm discussed in the previous section to separate the ID ((x) 
into Cl(x) and Cr^); which in turn separate tp(r) as (9(x,y,z) needs not to be separated): 
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|^(r)| 2 = %|^|^(r)| 2 
\Mr)\ 2 = l j$fmr)\* 

(21) 

Thus, in practice, if we already have the three dimensional eigen states {e, i/j}, all we need 
is to get the corresponding ID effective potential E(x). Note the current in x direction will 
be the same as from the ID formula because 8 2 (x,y, z) normalize to 1 over y and z at any 
x. We will ignore any current in the y and z directions. For the effective mass Hamiltonian, 
E(x) can be obtained from Eq.(17) as: 



' jmm* =3?w+*<*> <*) 

C. LCBB calculation for the 3D model 

So far, we have only used effective mass Hamiltonian to derive our TBS algorithm from 
ipi(r) to ipf(r) and ipf(r). However, we will like to use atomistic Hamiltonian (e.g.,EPM) to 
obtain ipi{r). One might ask whether the use of effective mass model derived formula will 
devalue the atomistic result into effective mass level. The answer is no. The reason is that 
the atomic features are built in ipi(r) by the atomistic Hamiltonian, while the separation 
algorithm is built on the smooth part (envelop function part) of ipi(r), which can be described 
by the effective mass formalism. Using the linear combination of bulk band method (LCBB 
method)^- to solve the electronic eigen states {Ei, ipi(r)}, the wave function tp{r) will be 
expanded by the Bloch states of the constituent bulk solids: 

^r) = Y J Cn,kU n , k {r)e tk - r (23) 

n,k 

where the periodic part u n ^{r) of the Bloch function is described by the plane wave 
functions as 

Un ^r) = -^=j2 A ^ G y G ' r ( 24 ) 
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Here, n is the band index, and k is the supercell reciprocal-lattice vector defined within 
the first BZ of the silicon primary cell while G is the reciprocal lattice of the primary cell 
chosen within an energy cutoff. The Hamiltonian matrix elements are evaluated within the 
basis set {u n ^(r)e lk ' r } , and the resulting Hamiltonian matrix is diagonalized to yield {C nt k}- 

We now go back to yield the ID potential E(x) in Eq(22) using the LCBB calculations, 
keep in mind that the separation formalisms in the subsections A and B correspond to the 
envelop function properties with the atomic features (e.g., u n ^{r)) of the wave functions and 
the potentials already being removed (or smoothed out). The one-dimensional integrated 
potential E v (x) in Eq.(22) can be calculated directly from the electrostatic potential <p{r) 
and the quantum confinement potential V str in Eq(l) as 

h v(x) = , (25) 

J \ip{r)\ 2 dydz 

The electrostatic potential caused by the carried charge density <f>(r) is generally smooth. 
The one-dimensional kinetic energy E v ^{x) is the effective mass kinetic energy, which involves 
the Laplacian on the envelope function of the wave function (this is different from the true 
kinetic anergy of the whole wave function as in Eq.(l)). In our LCBB representation, 
envelope function part of the wave function corresponds to the k component in Eq.(23), not 
the G component in Eq.(24). Thus if we define $(r) = (-^£2 ~ ^£2) \enveio P *P(r), then 
under LCBB expansion, we have: 

= E(^(^ - ^)) 2 + 2^(*« - Uk)?)C n , k u n ^ k - r (26) 

n,k y 2 

here ko(k) is the bottom of valley point for each k point. For example, in the Si case 
studied in the current paper, ko is the 0.83X point of the bottom of conduction band. Note 
that evaluating ip(r) adds no computational expense to the LCBB routine. Once we have 
■0( r )> the one-dimensional kinetic energy E y ^(x) can be calculated by 

_ fi>*(r)$(r)dydz 

J \^{ r )\Hydz [ n 

This concludes our full TBS algorithm using LCBB: first to calculate the eigen states 

{€i,ipi(r)} under LCBB, then use Eqs.(25)-(27) to obtain the ID potential E(x), then use 

Eqs.(20),(21) and ID formalism to separare ipi(r) into ^(r) and ^(r). The calculation 
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of the carrier charge density follow naturally from analogues equation of Eq(15), and the 
selfconsistent calculation is done via Eq.(2). Finally, the current is evaluated from Eq.(15). 



III. SIMULATION FOR 22-NM DOUBLE GATE ULTRA-THIN-BODY FIELD- 
EFFECT TRANSISTOR 

Following the ITRS 22 nm technology node, we present a three-dimensional atomistic 
quantum mechanical simulation on a 22 nm double gate ultra-thin-body field-effect transistor 
using the above presented three-dimensional TBS model. Previously, this kind of nanodevice 
has been studied by a TB and NEGF approach 41 . Here, to achieve a comparison between our 
TBS model and NEGF model, we use the same structural and material parameters as in Ref. 
41. Fig. 6 shows the structure and parameters of the simulated device. The channel direction 
is along (100). The gate length and silicon body thickness is Lq = 22 nm and t$i = 4.9 nm 
respectively. The thickness of the oxide is tox = 1.3 nm. The channel region is undoped 
while the source and drain region is highly n+ doped as Nd = 1 x 10 20 cm~ 3 . The drain bias 
is fixed to be Vd = 1-0V in the simulation. The Schrodinger equation (1) is solved in the 
whole device region with 0.8 million atoms using the LCBB method. An artificial periodic 
boundary condition is used connecting the left and right ends of the device with an artificial 
potential barrier between them. The resulting eigen states {E iy ipi{r)} are occupied by the 
source/drain Fermi- Levels by the three-dimensional algorithm in the previous section. Then 
the occupied electron density is used in the Poisson equation (2) to form a self-consistent 
calculation. The Poisson equation (2) is only solved within the box of the blue dashed line 
shown in Fig.6. A typical Pulay DIIS charge mixing iteration scheme is used to accelerate 
the speed of the self-consistent calculation. For more details for how to solve the Poisson 



equation and its boundary conditions, please see Refj34j. In the LCBB calculation, k-points 
are chosen from six X-valleys (two X100, two X010 and two X001) and in total there are 
about 7000 k-points included in the LCBB expansion. Two conduction bands are selected 
in the LCBB basis set. The calculation is carried out on a computer with a single CPU. 

Fig. 7 and 8 shows how the TBS works for such LCBB calculated million atom 3D 
system. For ballistic and tunneling case, the black solid line in Figs. 7 and 8 indicates the 
total wave of the eigen state while the red and blue lines with symbols indicate the right 
and left running (tunneling) waves. As can be seen from the figure, the two separated parts 
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coincide very well with the total wave function at both side of the barrier. For the ballistic 
case as shown in Fig. 7, the right running and left running waves are separated at the 
top of the barrier and then are injected ballistically into the other sides of the barrier. For 
the tunneling; shown in Fig. 8, the right and left tunneling waves separate at the 

minimum wave function point, then decay exponentially to the other side of the barrier. 
Note that unlike in the ID case shown in Figs.l to 3, the potential profile E(x) in the 3D 
case, hence the potential maximum point x m depends on the state index i. Nevertheless, 
the separating algorithm remains the same. 

Fig. 9 shows the converged occupied electron density and local conduction band profile 
(which equals 0(r) + -E^V) f° r Vg = 0.0 — 0.8V. As can be seen from the figure, in low gate 
voltage cases, the electron density decays very fast from the source/drain to the channel. 
The electron density at small gate bias decays 10 7 order from the source/drain region to 
the channel region. As the gate bias gets higher, the barrier in the channel is pushed down 
leading to a significantly increased charge density in the channel. It should be noted that the 
maximum point of the barrier moves towards the source side as the gate potential further 
pushes the barrier down. 

Finally, Fig. 10 compares the I/V curves of our splitting model and the NEGF incor- 
porated sp 3 d 5 s* tight-binding model^^i. The I/V data of the TB+NEGF model is from 
Ref.41. Despite of the different Hamiltonians used, and different numerical details in solving 
the Poisson equations, the results aggree excellently. Tabled gives a detailed comparison 
of some key device performance parameters. As can be seen from the table, there is just 
a 12.7% difference between the ON currents of the two models and only 10 mV differ- 
ence between the threshold voltages. Both two models give exactly the same sub-threshold 
swing S = 63mV/dec. which is defined as S = dL ^" oId i n sub-threshold region. This value 
(S = 63mV/dec.) is close to the theoretical limit of kT ■ InlO = 60meV, indicating a very 
effective gate control in such a double-gate ultra-thin-body device. 

Note that the TB+NEGF model solves the transport problem based on an open boundary 
condition while our splitting model is based on a periodic boundary condition and solves the 
same problem from the eigen states. From the comparisons, we see that these two methods 
lead to almost the same results. However, the computational costs for these two models are 
significantly different, while thousands of computer processors are used for the TB+NEGF 
approach, a single processor is used for our TBS approach. Overall, this is an evidence of 
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the validity of our TBS model for 3D systems. 
IV. SUMMARY 

In summary, we have presented a new empirical pseudopotential calculation approach 
to simulate million atom nanostructure semiconductor devices using the periodic boundary 
conditions. To treat the non-equilibrium condition, instead of calculating the scattering 
states from the source and drain, we calculated the stationary states by the linear combina- 
tion of bulk band method and then separated the whole wave function into source and drain 
parts according to a top of the barrier splitting (TBS) scheme based on the physical insight 
of ballistic and tunneling transport. The separated electronic states were then occupied 
according to the source/drain Fermi-Levels to yield the occupied electron density which is 
then used in the Poisson Equation solver to form a self-consistent calculation. The validity 
of TBS was verified by a comparison between the TBS calculation and scattering states 
calculation in a ID effective mass model. It is also verified by comparing our LCBB-TBS 
calculation for a 3D double-gate field-effect transistor with a non-equilibrium Green's func- 
tion tight-binding result. However, instead of using thousands of processors, our method can 
be calculated using a single processor. Thus, this method provides a fast way to simulate 
the non-equilibrium million atom nanostructure problems. 
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FIG. 1: The potential barrier and different energy levels. 





TBS 


TB+NEGF 


Ion (fiA/fim) 


3265 


3740 


V th (mV) 


460 


450 


S (mV/dec.) 


63 


63 



TABLE I: Comparison of some key parameters of device performance with Ref.41. 
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FIG. 2: The splitted wave functions for the ballistic case (a) and compared with the direct 
calculated scattering states (b). 
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FIG. 3: The splitted wave functions for the tunneling case and compared with the direct calculated 
scattering states. 



23 



T 




x position (a.u.) 

FIG. 4: The occupied charge density using equation (14) with the splitted wave functions in com- 
parison with the charge density calculated with directly calculated scattering state using equation 
(3). The oscillation at the two ends is due to the artificial periodic boundary condition used to 
calculate ipi(x). 
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FIG. 5: The total current as a function of gate potential V g as calculated from equation (14) in 
comparison with the current calculated from equation (4) using scattering states. The ft is an 
artificially sharp truncation for the occupation function to test the tunneling current. 
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FIG. 6: (Color online) Structure of the simulated 22 nm Double Gate Ultra-Thin-Body Filed-Effect 
Transistor. 
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FIG. 7: (Color online) The splitted wave functions for ballistic case. 
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FIG. 8: (Color online) The splitted wave functions for tunneling case. 
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FIG. 9: (Color online) ID profile of the occupied electron density (a) and the conduction band 
profile (b) for gate bias V g = 0.0 - 0.81/. 




FIG. 10: (Color online) Transfer characteristics compared with Ref.41 from the 3D quantum 
mechanical model. 
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